#####################################################################
##Replication Code for ##############################################
##The Politics of Promotion in China's Foreign Policy Bureaucracy ###
##Tyler Jost and Yucong Li ##########################################
##The China Quarterly ###############################################
##(4) Appendix Analyses #############################################
#####################################################################

#########################################
####Load Data############################
#########################################
rm(list = ls())
load("ts_minister.RData")
load("ts_minister_connections.RData")
load("ts_vminister.RData")
load("ts_vminister_connections.RData")
load("mfa_appointments.RData")
load("mfa_personnel.RData")
load("descriptives.RData")

crse <- function(x, df){
  d.cse <- na.omit(df[, c("id", all.vars(formula(x)))])
  fc <- d.cse$id
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <- coeftest(x, c.vcov)
  return(result)
}

crse2 <- function(x, df){
  d.cse <- na.omit(df[, c("year", all.vars(formula(x)))])
  fc <- d.cse$year
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <- coeftest(x, c.vcov)
  return(result)
}

##2 Additional Analyses#######################
##############################################

##2.2 Cross-Sectional Analysis################
##############################################
##Cross-Sectional of Ambassadorial Pool#######
##############################################

cross <- vmin_connections %>%
  arrange(cname, desc(year))

cross <- cross %>%
  distinct(cname, .keep_all = TRUE)

# Merge the ildstart column from vmin into cross based on cname
cross <- cross %>%
  mutate(ildstart = vmin$ildstart[match(cname, vmin$cname)])

## Table 3
m1 <- glm(seniorshareabroad ~ juniorhomeposts + juniorabroadposts + juniorpriorityposts + as.factor(pc) + as.factor(count_bin), 
          data=cross)
crse(m1, cross)
cluster1 <- crse(m1, cross)
-0.0283055*sd(cross$juniorhomeposts,na.rm = T) #8% decrease
sd(cross$juniorhomeposts,na.rm = T)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(seniorshareabroad ~ juniorhomeposts + juniorabroadposts + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=cross)
crse(m2, cross)
cluster2 <- crse(m2, cross)
pseudoR2 <- round(pR2(m2)[[4]],2)
-0.02514901*sd(cross$juniorhomeposts,na.rm = T) #7% decrease

stargazer(m1, m2, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Cross-Sectional of Ambassadorial Pool"),
          covariate.labels = c("Junior Home Posts (count)", "Junior Abroad Posts (count)", "Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2)),
          dep.var.labels = c("Senior Share of Time Abroad"),
          title = "Cross-Sectional Analysis of Ambassadorial Pool",
          notes.append=FALSE, no.space=T,
          notes = c("Heteroskedasticity-Robust standard errors.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

##############################################
##Promotion Analysis##########################
##############################################
##Cross-Sectional of Ambassadorial Pool#######
##############################################

## Table 4
m1 <- glm(vm_promote ~ seniorshareabroad + as.factor(pc) + as.factor(count_bin), data=cross, family="binomial")
crse(m1, cross)
cluster1 <- crse(m1, cross)
pseudoR1 <- round(pR2(m1)[[4]],2)
beta <- coef(m1)["seniorshareabroad"]
rounded_percentage_change <- round((exp(beta * sd(cross$seniorshareabroad)) - 1) * 100, 2)
rounded_percentage_change

m2 <- glm(vm_promote ~ seniorshareabroad  + midcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
              data=cross, family="binomial")
crse(m2, cross)
cluster2 <- crse(m2, cross)
pseudoR2 <- round(pR2(m2)[[4]],2)

stargazer(m1, m2, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Cross-Sectional of Ambassadorial Pool"),
          covariate.labels = c("Senior Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)","Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2)),
          dep.var.labels = c("Promotion to Vice Minister Rank"),
          title = "Cross-Sectional Analysis of Ambassadorial Pool",
          notes.append=FALSE, no.space=T,
          notes = c("Heteroskedasticity-Robust standard errors.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))



##############################################
##2.3 Interaction: MFA Professionalism########
##############################################
d.stat <- descriptives %>% 
  group_by(year) %>% 
  summarise(male = mean(male, na.rm = T),
            long_march = mean(long_march, na.rm = T),
            military = mean(military, na.rm = T),
            civilcollege = mean(civilcollege, na.rm = T),
            abroad = mean(abroad, na.rm = T),
            ild = mean(ild, na.rm = T))
d.stat <- as.data.frame(d.stat)
professional <- subset(d.stat, select=c("year", "military"))
colnames(professional) <- c("year", "sharemilitary")

vmin <- merge(vmin, professional, by="year")
#Rescale for interpretability
vmin$sharemilitary <- 1 - vmin$sharemilitary

## Table 5
interaction <- glm(vm_promote ~ seniorshareabroad*sharemilitary + midcount + totaltreaty + juniorpriorityposts + male + military +ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
                   data=vmin, family="binomial")
crse(interaction, vmin)
cluster1 <- crse(interaction, vmin)
pseudoR1 <- round(pR2(interaction)[[4]],2)

stargazer(interaction, type = "latex", 
          se = list(cluster1[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Interaction Term between Senior Share Abroad and Share without Military Background"),
          covariate.labels = c("Senior Share of Time Abroad", "Share of MFA without Military Background", "International Disputes (count)", "Diplomatic Treaties (count)",  "Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD" , "Higher Civilian Education", "Princeling", "Senior Share of Time Abroad*Share without Military Background"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1)),
          dep.var.labels = c("Promotion to Vice Minister Rank"),
          title = "Effect of MFA without Military Background and Senior Share of Time Abroad on Vice Minister Promotions",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

## Figure 2
p1 <- interplot(interaction, "seniorshareabroad", "sharemilitary") + 
  theme_minimal() + xlab("Share of MFA without Military Background") + ylab("Effect of Higher Share of Time Abroad") + scale_y_continuous(limits = c(-4, 0)) + scale_x_continuous(limits = c(0.2, 1)) +
  theme(legend.position="bottom", legend.margin=margin(-12.5,0,0,0), legend.text = element_text(face="bold"),
        axis.title.y = element_text(face="bold"), legend.title=element_blank(), 
        axis.text=element_text(face="bold"), text = element_text(size = 10))

ggsave("professionalism_posting_interaction.pdf", plot = p1, scale = 1, width = 4.5, height = 3, dpi = 500, limitsize = T)


##############################################
##2.4 Interaction: Domestic Instability#######
##############################################

unstable_years <- c(1958:1962, 1966:1976, 1989, 2012:2013)
vmin$domestic_instability <- ifelse(vmin$year %in% unstable_years,1,0)

## Table 6
interaction <- glm(vm_promote ~ seniorshareabroad*domestic_instability + midcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
                    data=vmin, family="binomial")
crse(interaction, vmin)
cluster1 <- crse(interaction, vmin)
pseudoR1 <- round(pR2(interaction)[[4]],2)

stargazer(interaction, type = "latex", 
          se = list(cluster1[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Interaction Term between Foreign Postings and Domestic Stability"),
          covariate.labels = c("Senior Share of Time Abroad", "Domestic Instability", "International Disputes (count)", "Diplomatic Treaties (count)", "Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling","Senior Share of Time Abroad*Domestic Instability "),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1)),
          dep.var.labels = c("Promotion to Vice Minister Rank"),
          title = "Impact of Foreign Postings and Domestic Stability on Vice Minister Promotions",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))


## Figure 3
p1 <- interplot(interaction, "seniorshareabroad", "domestic_instability") + 
  theme_minimal() + xlab("High Domestic Instability") + ylab("Effect of Higher Share of Time Abroad") + scale_y_continuous(limits = c(-4, 1)) +
  theme(legend.position="bottom", legend.margin=margin(-12.5,0,0,0), legend.text = element_text(face="bold"),
        axis.title.y = element_text(face="bold"), legend.title=element_blank(), 
        axis.text=element_text(face="bold"), text = element_text(size = 10))
ggsave("stability_posting_interaction.pdf", plot = p1, scale = 1, width = 4.5, height = 3, dpi = 500, limitsize = T)

##############################################
##2.5 Social Ties Analysis####################
##############################################
##Cross-Sectional of Ambassadorial Pool#######
##############################################

## Table 7
m1 <- glm(totalconnections_count ~ juniorhomeposts + juniorabroadposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
            data=cross)
crse(m1, cross)
cluster1 <- crse(m1, cross)
pseudoR1 <- round(pR2(m1)[[4]],2)
3.90428*sd(cross$juniorhomeposts,na.rm = T) 
2.20841*sd(cross$juniorabroadposts,na.rm = T) 


m2 <- glm(totalconnections_count ~ seniorshareabroad + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
            data=cross)
crse(m2, cross)
cluster2 <- crse(m2, cross)
pseudoR2 <- round(pR2(m2)[[4]],2)
-23.29734*sd(cross$seniorshareabroad, na.rm = T) 
sd(cross$seniorshareabroad)

stargazer(m1, m2, type = "latex", 
          se = list(cluster1[,"Std. Error"],cluster2[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Cross-Sectional Analysis of Social Ties"),
          covariate.labels = c("Junior Home Posts (count)", "Junior Abroad Posts (count)", "Senior Share of Time Abroad", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2)),
          dep.var.labels = c("Total Connections"),
          title = "Cross-Sectional Analysis of Social Ties",
          notes.append=FALSE, no.space=T,
          notes = c("Heteroskedasticity-Robust standard errors.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

######################################################################
##2.6 Alternative measurement strategies of international disputes####
######################################################################

## Table 8
##Vice Minister Level Promotion#########################
m1 <- glm(vm_promote ~ seniorshareabroad + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster1 <- crse(m1, vmin)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidsuccesscount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster2 <- crse(m2, vmin)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ seniorshareabroad + totalmidsuccesscount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster3 <- crse(m3, vmin)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ seniorshareabroad + totalmidsuccesscount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster4 <- crse(m4, vmin)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ seniorshareabroad + totalmidsuccesscount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster5 <- crse(m5, vmin[vmin$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ seniorshareabroad + totalmidsuccesscount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster6 <- crse(m6, vmin[vmin$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad", "China-Favored Dispute Outcomes (count)", "Diplomatic Treaties (count)",
                               "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

##Minister Level Promotion#########################
m1 <- glm(high_promote ~ seniorshareabroad + as.factor(count_bin), 
          data=minister, family="binomial")
cluster1 <- crse2(m1, minister)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ totalmidsuccesscount + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster2 <- crse2(m2, minister)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ seniorshareabroad + totalmidsuccesscount + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster3 <- crse2(m3, minister)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ seniorshareabroad + totalmidsuccesscount + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister, family="binomial")
cluster4 <- crse2(m4, minister)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ seniorshareabroad + totalmidsuccesscount + alltreaties + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster5 <- crse2(m5, minister[minister$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ seniorshareabroad + totalmidsuccesscount + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster6 <- crse2(m6, minister[minister$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

## No one promoted to Minister Level had successful mid experience
minister %>%
  dplyr::select(cname, high_promote, totalmidsuccesscount) %>%
  distinct() %>%
  count(high_promote, totalmidsuccesscount) %>%
  print(n = Inf)

###############################################################
## 2.7 Alternative Measurement Strategies of Treaties##########
#######Treaty with Major Power (Minister Promotion Only) ######
###############################################################
m1 <- glm(high_promote ~ seniorshareabroad + as.factor(count_bin), 
          data=minister, family="binomial")
cluster1 <- crse2(m1, minister)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ allmid + totalprioritytreaty + as.factor(count_bin), 
          data=minister, family="binomial")
cluster2 <- crse2(m2, minister)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ seniorshareabroad + allmid + totalprioritytreaty + as.factor(count_bin), 
          data=minister, family="binomial")
cluster3 <- crse2(m3, minister)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ seniorshareabroad + allmid + totalprioritytreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister, family="binomial")
cluster4 <- crse2(m4, minister)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ seniorshareabroad + allmid + totalprioritytreaty + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster5 <- crse2(m5, minister[minister$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ seniorshareabroad + allmid + totalprioritytreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster6 <- crse2(m6, minister[minister$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)


## No one promoted to minister level had signed treaties with major power in senior experience 
minister %>%
  dplyr::select(cname, high_promote, totalprioritytreaty) %>%
  distinct() %>%
  filter(high_promote == 1) %>%
  count(totalprioritytreaty)%>%
  print(n = Inf)

##2.8 Spokesperson Analysis ########################
####################################################
##Difference in Means of Spokesperson###############
####################################################
spokesperson.names <- c( "钱其琛","齐怀远","俞志忠","王振宇", "马毓真", "李肇星","李金华","金桂华","段津","吴建民","范慧娟","李建英","陈健","沈国放","崔天凯","唐国强","朱邦造",
                         "孙玉玺", "章启月", "孔泉", "刘建超", "秦刚", "姜瑜", "马朝旭", "洪磊", "刘为民", "华春莹", "陆慷", "耿爽", "赵立坚", "汪文斌", "毛宁", "林剑")
cross <- cross %>%
  mutate(spokesperson = ifelse(cname %in% spokesperson.names, 1, 0))

## Separate individuals who held the spokesperson position and those that did not
group_spokesperson <- cross %>% 
  filter(spokesperson == 1, year >= 1982) %>% 
  pull(vm_promote)

group_non_spokesperson <- cross %>% 
  filter(spokesperson == 0, year >= 1982) %>% 
  pull(vm_promote)


## Conduct a t-test
t_test_result <- t.test(group_spokesperson, group_non_spokesperson, var.equal = TRUE)
t_test_result

## Visualization 
promotion_summary <- cross %>%
  group_by(spokesperson) %>%
  summarize(
    promotion_rate = mean(vm_promote, na.rm = TRUE),
    count = n(),
    se = sd(vm_promote, na.rm = TRUE) / sqrt(count) 
  )

promotion_summary <- promotion_summary %>%
  mutate(
    ci_lower = promotion_rate - qt(0.975, count - 1) * se,
    ci_upper = promotion_rate + qt(0.975, count - 1) * se
  )


## Figure 4
barplot <- ggplot(promotion_summary, aes(x = as.factor(spokesperson), y = promotion_rate)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  labs(x = "Spokesperson Experience (0 = No, 1 = Yes)", 
       y = "Mean Promotion Rate") +
  theme_minimal() +
  theme(
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    plot.title = element_text(size = 14, face = "bold")
  )

ggsave("spokesperson.pdf", plot = barplot, scale = 1, width = 10, height = 4, dpi = 500, limitsize = T)

###############################################################
## Cross-sectional model for spokesperson #####################
###############################################################

## Table 9
m1 <- glm(vm_promote ~ seniorshareabroad  + spokesperson + midcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
            data=cross[cross$year>=1982,], family="binomial")
robust_se_1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC0"))
pseudoR1 <- round(pR2(m1)[[4]],2)

stargazer(m1, type = "latex", 
          se = list(robust_se_1[,"Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Cross-Sectional of Ambassadorial Pool"),
          covariate.labels = c("Senior Share of Time Abroad", "MFA Spokesperson","International Disputes (count)", "Diplomatic Treaties (count)","Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR1)),
          dep.var.labels = c("Promotion to Vice Minister Rank"),
          title = "Cross-Sectional Analysis of Ambassadorial Pool",
          notes.append=FALSE, no.space=T,
          notes = c("Heteroskedasticity-Robust standard errors.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

###############################################################
## Interaction between Spokesperson and MIDs ##################
###############################################################
mids <- read.csv("dyadic_mid_4.02.csv")
mids <- mids[mids$statea==710 & mids$year>=1949,]
mids <- subset(mids, select=c("stateb", "endyear", "outcome"))
mids$favorablemidout <- ifelse(mids$outcome %in% c(1,4,6),1,0)
counts <- mids %>%
  group_by(endyear) %>%
  summarize(
    midcount = n(),
    midsuccesscount = sum(favorablemidout == 1, na.rm = TRUE)
  ) %>%
  ungroup()
counts <- as.data.frame(counts)
colnames(counts) <- c("year", "midtotal", "midsuccesstotal")

## Filter and expand each spokesperson's term to individual years

spokesperson_years <- d %>%
  filter(clean_title == "MFA Spokesperson") %>%
  rowwise() %>%
  mutate(year = list(seq(start, end))) %>%
  unnest(cols = year) 
spokesperson_years$incumb <- 1
spokesperson_years <- spokesperson_years[, c("cname", "year", "incumb")]

## Count all MIDs involving China each year for each spokesperson
spokesperson_years <- spokesperson_years %>%
  left_join(counts, by = "year") 

## Replace NA values in midtotal and midsuccesstotal with 0
spokesperson_years <- spokesperson_years %>%
  mutate(
    midtotal = replace_na(midtotal, 0),
    midsuccesstotal = replace_na(midsuccesstotal, 0)
  )


## Overall disputes handled while in the position of spokesperson
spokesperson_years <- spokesperson_years %>%
  group_by(cname) %>%
  arrange(cname, year) %>%
  mutate(
    midtotalsum = sum(midtotal),
    midsuccesstotalsum = sum(midsuccesstotal)
  ) %>%
  ungroup()

spokesperson_years <- spokesperson_years %>%
  distinct(cname, .keep_all = TRUE)


## Merge the overall number of disputes into cross sectional df
cross <- cross %>%
  left_join(spokesperson_years, by = "cname")

## Replace NA values in midtotal and midsuccesstotal with 0
cross <- cross %>%
  mutate(
    midtotalsum = replace_na(midtotalsum, 0),
    midsuccesstotalsum = replace_na(midsuccesstotalsum, 0)
  )

## Update midtotalsum and midsuccesstotalsum when spokesperson == 0 (ambassador only has bilatery while spokesperson has exposure to all)
cross <- cross %>%
  mutate(
    midtotalsum = ifelse(spokesperson == 0, totalmidcount, midtotalsum),
    midsuccesstotalsum = ifelse(spokesperson == 0, totalmidsuccesscount, midsuccesstotalsum)
  )


## Interaction Model with cross-sectional df

## Table 10
interaction2 <- glm(vm_promote ~ seniorshareabroad + spokesperson*midtotalsum + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
                    data=cross[cross$year.x>=1982,], family="binomial")
interaction2
robust_se_2 <- coeftest(interaction2, vcov = vcovHC(interaction2, type = "HC0"))
pseudoR2 <- round(pR2(interaction2)[[4]],2)

interaction3 <- glm(vm_promote ~ seniorshareabroad + spokesperson*midsuccesstotalsum + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
                    data=cross[cross$year.x>=1982,], family="binomial")
interaction3
robust_se_3 <- coeftest(interaction3, vcov = vcovHC(interaction3, type = "HC0"))
pseudoR3 <- round(pR2(interaction3)[[4]],2)

stargazer(interaction2, interaction3, type = "latex", 
          se = list(robust_se_2[, "Std. Error"], robust_se_3[, "Std. Error"]),
          omit=c("Constant","pc","count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: Cross-Sectional of Ambassadorial Pool"),
          covariate.labels = c("Senior Share of Time Abroad", "Spokesperson", "International Disputes (count)", "China-Favored Dispute Outcomes (count)", "Diplomatic Treaties (count)","Junior Priority Posts (count)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$"),
                         c("McFadden Pseudo-R$^2$", pseudoR2, pseudoR3)),
          dep.var.labels = c("Promotion to Vice Minister Rank"),
          title = "Cross-Sectional Analysis of Ambassadorial Pool",
          notes.append=FALSE, no.space=T,
          notes = c("Heteroskedasticity-Robust standard errors.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

## Figure 5
p2 <- interplot(interaction2, "midtotalsum", "spokesperson") + 
  theme_minimal() + xlab("Holding Spokesperson Position (0 = No, 1 = Yes)") + ylab("International Disputes (count)") + 
  theme(legend.position="bottom", legend.text = element_text(face="bold"),
        axis.title.y = element_text(face="bold"), legend.title=element_blank(), 
        axis.text=element_text(face="bold"), text = element_text(size = 10))

ggsave("spokespersontreaty1.pdf", plot = p2, scale = 1, width = 4.5, height = 3, dpi = 500, limitsize = T)

## Figure 6
p3 <- interplot(interaction3, "midsuccesstotalsum", "spokesperson") + 
  theme_minimal() + xlab("Holding Spokesperson Position (0 = No, 1 = Yes)") + ylab("China-Favored Dispute Outcomes (count)") + 
  theme(legend.position="bottom", legend.text = element_text(face="bold"),
        axis.title.y = element_text(face="bold"), legend.title=element_blank(), 
        axis.text=element_text(face="bold"), text = element_text(size = 10))

ggsave("spokespersontreaty2.pdf", plot = p3, scale = 1, width = 4.5, height = 3.5, dpi = 500, limitsize = T)


###############################################################
##2.9 Major IOs experience and general abroad experience ######
###############################################################

## UN Experiences
vmin <- vmin %>%
  group_by(cname) %>%
  mutate(ioexperience = cumsum(expertise == "Multilateral Organization" & location != "APEC"),
         nonioabroadshare = (seniorabroadexperience - ioexperience) / (seniorabroadexperience + beijingexperience),
         ioexperienceshare = ioexperience / ((seniorabroadexperience + beijingexperience))) %>%
  ungroup()


minister <- minister %>%
  group_by(cname) %>%
  mutate(ioexperience = cumsum(expertise == "Multilateral Organization" & location != "APEC"),
         nonioabroadshare = (seniorabroadexperience - ioexperience) / (seniorabroadexperience + beijingexperience),
         ioexperienceshare = ioexperience / ((seniorabroadexperience + beijingexperience))) %>%
  ungroup()

#############################
## General IOs experience####
#############################

IOs <- c("APEC", "ASEAN", "ECOWAS", "European Union ", "UN - Geneva", "UN - NYC", "UN - Vienna")

vmin <- vmin %>%
  group_by(cname) %>%
  mutate(io2experience = cumsum(location %in% IOs),
         nonio2abroadshare = (seniorabroadexperience - io2experience) / (seniorabroadexperience + beijingexperience),
         io2experienceshare = io2experience / ((seniorabroadexperience + beijingexperience))) %>%
  ungroup()

minister <- minister %>%
  group_by(cname) %>%
  mutate(io2experience = cumsum(location %in% IOs),
         nonio2abroadshare = (seniorabroadexperience - io2experience) / (seniorabroadexperience + beijingexperience),
         io2experienceshare = io2experience / ((seniorabroadexperience + beijingexperience))) %>%
  ungroup()

## Table 11
## VM Level UN experience
m1 <- glm(vm_promote ~ ioexperienceshare + nonioabroadshare + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster1 <- crse(m1, vmin)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster2 <- crse(m2, vmin)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ ioexperienceshare + nonioabroadshare + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster3 <- crse(m3, vmin)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ ioexperienceshare + nonioabroadshare + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster4 <- crse(m4, vmin)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ ioexperienceshare + nonioabroadshare + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster5 <- crse(m5, vmin[vmin$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ ioexperienceshare + nonioabroadshare + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster6 <- crse(m6, vmin[vmin$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad in UN", "Share of Time Abroad Non-UN", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

## Minister Level 
m1 <- glm(high_promote ~ ioexperienceshare + nonioabroadshare + as.factor(count_bin), 
          data=minister, family="binomial")
cluster1 <- crse2(m1, minister)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster2 <- crse2(m2, minister)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ ioexperienceshare + nonioabroadshare + allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster3 <- crse2(m3, minister)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ ioexperienceshare + nonioabroadshare + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister, family="binomial")
cluster4 <- crse2(m4, minister)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ ioexperienceshare + nonioabroadshare + allmid + alltreaties + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster5 <- crse2(m5, minister[minister$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ ioexperienceshare + nonioabroadshare + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster6 <- crse2(m6, minister[minister$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

## No one promoted to minister level has UN senior experience
minister %>%
  dplyr::select(cname, high_promote, ioexperienceshare) %>%
  distinct() %>%
  count(high_promote, ioexperienceshare)

## Table 12
## VM Level General IOs Experience
m1 <- glm(vm_promote ~ io2experienceshare + nonio2abroadshare + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster1 <- crse(m1, vmin)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster2 <- crse(m2, vmin)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ io2experienceshare + nonio2abroadshare + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster3 <- crse(m3, vmin)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ io2experienceshare + nonio2abroadshare + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster4 <- crse(m4, vmin)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ io2experienceshare + nonio2abroadshare + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster5 <- crse(m5, vmin[vmin$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ io2experienceshare + nonio2abroadshare + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster6 <- crse(m6, vmin[vmin$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad in IOs", "Share of Time Abroad in Embassy", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))


## Minister Level  IOs Experience
m1 <- glm(high_promote ~ io2experienceshare + nonio2abroadshare + as.factor(count_bin), 
          data=minister, family="binomial")
cluster1 <- crse2(m1, minister)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster2 <- crse2(m2, minister)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ io2experienceshare + nonio2abroadshare + allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster3 <- crse2(m3, minister)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ io2experienceshare + nonio2abroadshare + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister, family="binomial")
cluster4 <- crse2(m4, minister)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ io2experienceshare + nonio2abroadshare + allmid + alltreaties + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster5 <- crse2(m5, minister[minister$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ io2experienceshare + nonio2abroadshare + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster6 <- crse2(m6, minister[minister$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

## No one promoted to minister level has IO senior experience
minister %>%
  dplyr::select(cname, high_promote, io2experienceshare) %>%
  distinct() %>%
  count(high_promote, io2experienceshare)


###############################################################
##2.10 Dropping Active Diplomats ##############################
###############################################################

## Filter rows where end == 2023 and extract unique cnames
active_list <- d %>%
  filter(end == 2023) %>%
  pull(cname) %>%
  unique()
active_list

vm_list <- vmin %>%
  filter(vm_promote == 1) %>%
  pull(cname)%>%
  unique
vm_list

## Remove "秦刚" from the list (the only one who are not active in the list)
active_list <- setdiff(active_list, "秦刚")

## Check if names in vmin$cname exist in active_list who have not yet been promoted to vm level
vmin <- vmin %>%
  mutate(activediplomats = ifelse(cname %in% active_list & !(cname %in% vm_list), 1, 0))

## Create the new dataframe, keeping only ex diplomats
vminended <- vmin %>%
  filter(activediplomats == 0)

## Check if names in minister$cname exist in active_list who have not yet been promoted to minister level
min_list <- minister %>%
  filter(high_promote == 1) %>%
  pull(cname)%>%
  unique
min_list

minister <- minister %>%
  mutate(activediplomats = ifelse(cname %in% active_list & !(cname %in% min_list), 1, 0))

# minister <- minister %>%
#   mutate(activediplomats = ifelse(cname %in% active_list, 1, 0))

ministerended <- minister %>%
  filter(activediplomats == 0)

## Redo the promotion analysis

## Table 13
# VM promotion
m1 <- glm(vm_promote ~ seniorshareabroad + as.factor(pc) + as.factor(count_bin), 
          data=vminended, family="binomial")
cluster1 <- crse(m1, vminended)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vminended, family="binomial")
cluster2 <- crse(m2, vminended)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vminended, family="binomial")
cluster3 <- crse(m3, vminended)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vminended, family="binomial")
cluster4 <- crse(m4, vminended)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vminended[vminended$year>=1982,], family="binomial")
cluster5 <- crse(m5, vminended[vminended$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vminended[vminended$year>=1982,], family="binomial")
cluster6 <- crse(m6, vminended[vminended$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("The analysis excludes diplomats who are currently active in the Ministry of Foreign Affairs (MFA) but have not yet been promoted to the Vice Minister level.", 
                    "Robust standard errors are clustered by individual. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01"))

## Table 14
# Minister level promotion
m1 <- glm(high_promote ~ seniorshareabroad + as.factor(count_bin), 
          data=ministerended, family="binomial")
cluster1 <- crse2(m1, ministerended)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ allmid + alltreaties + as.factor(count_bin), 
          data=ministerended, family="binomial")
cluster2 <- crse2(m2, ministerended)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + as.factor(count_bin), 
          data=ministerended, family="binomial")
cluster3 <- crse2(m3, ministerended)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=ministerended, family="binomial")
cluster4 <- crse2(m4, ministerended)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + as.factor(count_bin), 
          data=ministerended[ministerended$year>=1982,], family="binomial")
cluster5 <- crse2(m5, ministerended[ministerended$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=ministerended[ministerended$year>=1982,], family="binomial")
cluster6 <- crse2(m6, ministerended[ministerended$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"],
                    cluster4[,"Std. Error"], cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: min_promotion"),
          covariate.labels = c("Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)", "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(minister$year)),3),
                           length(unique(minister$year[is.na(minister$juniorpriorityposts)==F])),
                           length(unique(minister$year[minister$year>=1982])),
                           length(unique(minister$year[is.na(minister$juniorpriorityposts)==F & minister$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Appointment to Minister Rank"),
          title = "Postings, Performance, and Appointment to Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Analysis excludes currently active MFA diplomats who have not yet been promoted to Minister level.", 
                    "Robust standard errors are clustered by year. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01"))

###############################################################
##2.11 FM Fixed Effect ########################################
###############################################################

vmin$fm <- NA
vmin$fm[vmin$year >= 1949 & vmin$year <= 1958] <- "Zhou Enlai"
vmin$fm[vmin$year > 1958 & vmin$year <= 1972] <- "Chen Yi"
vmin$fm[vmin$year > 1972 & vmin$year <= 1974] <- "Ji Pengfei"
vmin$fm[vmin$year > 1974 & vmin$year <= 1976] <- "Qiao Guanhua"
vmin$fm[vmin$year > 1976 & vmin$year <= 1982] <- "Huang Hua"
vmin$fm[vmin$year > 1982 & vmin$year <= 1988] <- "Wu Xueqian"
vmin$fm[vmin$year > 1988 & vmin$year <= 1998] <- "Qian Qichen"
vmin$fm[vmin$year > 1998 & vmin$year <= 2003] <- "Tang Jiaxuan"
vmin$fm[vmin$year > 2003 & vmin$year <= 2007] <- "Li Zhaoxing"
vmin$fm[vmin$year > 2007 & vmin$year <= 2013] <- "Yang Jiechi"
vmin$fm[vmin$year > 2013 & vmin$year <= 2022] <- "Wang Yi"
vmin$fm[vmin$year > 2022 & vmin$year <= 2023] <- "Qin Gang"

## Table 15
## Rerun VM promotion analysis including FM fixed effect
m1 <- glm(vm_promote ~ seniorshareabroad + as.factor(fm) + as.factor(count_bin), 
          data=vmin, family="binomial")
m1
cluster1 <- crse(m1, vmin)
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidcount + totaltreaty + as.factor(fm) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster2 <- crse(m2, vmin)
pseudoR2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(fm) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster3 <- crse(m3, vmin)
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(fm) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster4 <- crse(m4, vmin)
pseudoR4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(fm) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster5 <- crse(m5, vmin[vmin$year>=1982,])
pseudoR5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(fm) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster6 <- crse(m6, vmin[vmin$year>=1982,])
pseudoR6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "fm", "count_bin"),
          omit.stat=c("f", "ser"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Priority Experience (junior)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Foreign Minister Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4, pseudoR5, pseudoR6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))


####################################################################
##2.12 Additional Model with Current Beijing Assignment ############
####################################################################

## Table 16
m1 <- glm(vm_promote ~ beijing + seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year >= 1949 & vmin$year <= 1976, ], family="binomial")
cluster1 <- crse(m1, vmin[vmin$year >= 1949 & vmin$year <= 1976, ])
pseudoR1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ beijing + seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year >= 1949 & vmin$year <= 1976, ], family="binomial")
cluster2 <- crse(m2, vmin[vmin$year >= 1949 & vmin$year <= 1976, ])
pseudoR2 <- round(pR2(m2)[[4]],2)
vif(m2)

m3 <- glm(vm_promote ~ beijing + seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year >= 1976, ], family="binomial")
cluster3 <- crse(m3, vmin[vmin$year >= 1976, ])
pseudoR3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ beijing + seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year >= 1976, ], family="binomial")
cluster4 <- crse(m4, vmin[vmin$year >= 1976, ])
pseudoR4 <- round(pR2(m4)[[4]],2)
vif(m4)

stargazer(m1, m2, m3, m4, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser",  "aic", "ll"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Current Beijing Assignment", "Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Junior Priority Experience (number of posts)", "Male", "Military Background", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", length(unique(vmin$cname[vmin$year<1976])), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year<1976])),
                           length(unique(vmin$cname[vmin$year>=1976])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1976]))),
                         c("McFadden Pseudo-R$^2$", pseudoR1, pseudoR2, pseudoR3, pseudoR4)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

##############################################
##3 Diplomatic Removal Analysis###############
##############################################

##Any Position
purged <- individual[individual$purged==1 & is.na(individual$purged)==F,]
length(unique(purged$cname)) #some yrpurged are unavailable #118
purged$leader <- NA
purged$leader[purged$yrpurged<=1962] <- "Early Mao"
purged$leader[purged$yrpurged>1962 & purged$yrpurged<=1976] <- "Late Mao"
purged$leader[purged$yrpurged>1976 & purged$yrpurged<=1989] <- "Deng Xiaoping"
purged$leader[purged$yrpurged>1989 & purged$yrpurged<=2002] <- "Jiang Zemin"
purged$leader[purged$yrpurged>2002 & purged$yrpurged<=2012] <- "Hu Jintao"
purged$leader[purged$yrpurged>2012] <- "Xi Jinping"
purge1 <- data.frame(table(purged$leader))
colnames(purge1) <- c("Leader", "all")

##Assistant Minister/Ambassador + Above
seniorpurged <- individual[individual$seniorpurged==1 & is.na(individual$seniorpurged)==F,]
seniorpurged$leader <- NA
seniorpurged$leader[seniorpurged$yrseniorpurged<=1962] <- "Early Mao"
seniorpurged$leader[seniorpurged$yrseniorpurged>1962 & seniorpurged$yrseniorpurged<=1976] <- "Late Mao"
seniorpurged$leader[seniorpurged$yrseniorpurged>1976 & seniorpurged$yrseniorpurged<=1989] <- "Deng Xiaoping"
seniorpurged$leader[seniorpurged$yrseniorpurged>1989 & seniorpurged$yrseniorpurged<=2002] <- "Jiang Zemin"
seniorpurged$leader[seniorpurged$yrseniorpurged>2002 & seniorpurged$yrseniorpurged<=2012] <- "Hu Jintao"
seniorpurged$leader[seniorpurged$yrseniorpurged>2012] <- "Xi Jinping"
purge2 <- data.frame(table(seniorpurged$leader))
colnames(purge2) <- c("Leader", "Senior")
purge <- merge(purge1, purge2, by = "Leader", all.x = T, all.y = T)
purge$Junior <- purge$all - purge$Senior
purge <- subset(purge, select=c("Leader", "Senior", "Junior"))

table(seniorpurged$yrseniorpurged)

purge_long <- pivot_longer(purge, cols = c(Senior, Junior), names_to = "Position", values_to = "Count")
order <- c("Xi Jinping", "Hu Jintao", "Deng Xiaoping", "Late Mao", "Early Mao")

## Figure 7
barplot <- ggplot(purge_long, aes(y = factor(Leader, levels = order), x = Count, fill = Position, label = Count)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(position = position_dodge(0.9), hjust = -1, vjust = 0.5) +  # Adjust the position and alignment of labels
  labs(title = "", y = "", x = "") + scale_x_continuous(breaks = seq(0, 60, by = 10)) +
  theme_minimal() + scale_fill_manual(values = c("#FF5733", "#3399FF")) +
  theme(legend.position="bottom", legend.margin=margin(-6,0,0,0), legend.text = element_text(face="bold"),
        axis.title.y = element_text(face="bold"), legend.title=element_blank(), 
        axis.text=element_text(face="bold"), text = element_text(size = 18))

ggsave("removals.pdf", plot = barplot, scale = 1, width = 10, height = 4, dpi = 500, limitsize = T)

##Now, look at early depatures from position
d$duration <- d$end - d$start
earlyremove <- d[d$duration<2 & (d$clean_title=="Ambassador" | d$clean_title=="Assistant Minister" | d$clean_title=="Vice Minister") 
                            & d$exit_type!="Purge" & is.na(d$duration)==F,]
earlyremove <- earlyremove[is.na(earlyremove$cname)==F,]
length(unique(earlyremove$cname)) #160
length(unique(earlyremove$cname[earlyremove$exit_type=="Death"])) #9
length(unique(earlyremove$cname[earlyremove$exit_type=="Retire"])) #19
length(unique(earlyremove$cname[earlyremove$exit_type=="Transfer"])) #22
length(unique(earlyremove$cname[earlyremove$exit_type=="Within"])) #115
table(earlyremove$cname,earlyremove$exit_type)

